ctx

preamble

library(RANN)
library(arrow)
library(dplyr)
library(tidyr)
library(cluster)
library(circlize)
library(ComplexHeatmap)
library(SingleCellExperiment)
source("code/_utils.R")

loading

sce <- readRDS("outs/fil.rds")
ist <- readRDS("outs/lv1.rds")
jst <- readRDS("outs/lv2.rds")

labeling

sce$lv1 <- (kid <- ist$clust)[match(colnames(sce), names(kid))]
kid <- unlist(unname(lapply(jst, \(.) .$clust)))
idx <- match(colnames(sce), names(kid))
table(sce$lv2 <- kid[idx], exclude=NULL)

    BEC     epi     fib  FRCcts   FRCpv  FRCtcz     LEC    mCAF     AML      DC 
   1034     212     653     941     744    3753     945    2016    2943    1700 
 mac.pi  mono.c mono.nc     pDC     TAM     tum      NK  NK/ILC     Tcm     Tcn 
   1529    1276     213    1587     866   61871      42     422     197     747 
    Tex     Tfh     Tha     Thn      Tp    Treg 
   1969     921     310    3513     358    1071 

neighborhoods

cs <- split(colnames(sce), sce$sid)
names(id) <- id <- names(cs)
r <- 0.02; k <- 101
# get radial neighborhood
is <- lapply(id, \(sid) {
    sub <- sce[, cs[[sid]]]
    xy <- "Center(X|Y)_global_mm"
    xy <- grep(xy, names(colData(sub)))
    xy <- as.matrix(colData(sub)[xy])
    nn <- nn2(xy, searchtype="radius", r=r, k=k)
    is <- nn$nn.idx[, -1]; is[is == 0] <- NA; is
})

n <- unlist(lapply(is, \(.) rowSums(!is.na(.))))
l <- sprintf("#cells in %sum radial neighborhood", r*1e3)
hist(n, breaks=seq(-.5, max(n)+.5), xlab="", ylab="",  main=l)
abline(v=quantile(n, c(.05, .95)), col="red", lw=2, lty=2)
abline(v=m <- mean(n), col="red", lw=2); m

[1] 13.55114

composition

# quantify each neighborhood's
# subpopulation composition
fq <- lapply(id, \(.) {
    ks <- (se <- sce[, cs[[.]]])$lv2
    ks[is.na(ks)] <- "tum"
    mx <- matrix(ks[is[[.]]], nrow=ncol(se))
    mx[is.na(mx)] <- ""
    names(kt) <- kt <- unique(ks)
    ns <- sapply(kt, \(k) rowSums(mx == k))
    data.frame(
        prop.table(ns, 1),
        row.names=colnames(se))
}) |> bind_rows(); fq[is.na(fq)] <- 0

clustering

set.seed(251210)
km <- kmeans(fq, centers=k <- 10)$cluster
ns <- sprintf("%02d", seq_len(k))
ns <- factor(km, labels=paste0("N", ns))

visuals

sce$ctx <- ns[match(colnames(sce), names(ns))]
t(sapply(cs, \(.) table(ns[.]))); table(ns)

        N01  N02  N03 N04  N05  N06  N07 N08  N09 N10
011LNd    0 2545    0   0    1    0 1306 221  537  56
012LNd    0 4491    0   0    0    0  456   1   45   2
020LNm   13  956    0   0    2    0  929  72  319   4
031LNd    1 4214    0   0    1    0 2227  27  338   8
041SOd    0 2819    0   0    0    0  612  12  120   3
042SOd    0 3000    0   0    0    0  866  29  191   3
051LNn   88 1029    1   0   14    0 2694 314 2612 191
052LNn  147  460    4   0   12    0 1402 626 2486 674
061LNd    1 1437    0   0    1    0 1623 109  740  39
062LNd    0 2344    0   0    0    0 1110  14  198   1
071NPd  752    7 2314 213   12    0   15 117   62  60
072NPd 2287  397   70   0   19    0 1047 578 1502  41
080LNn    0 1892    0   0    2    0 2144 167  697 167
091EYn    0 1413    0   0    0    0  440  30   94   1
092EYn    0  222    0   0    4    0  324 227  338  86
100LNd    1 1122    0   0    4    0 1124 181  676  77
111LNn  114 1049    0   0   37    1 1204 814  906 195
112LNn   20 1110    0   0  110    0 1089 787  532  95
121LNm    5  306    0   0 2525 3435  589 486  282 330
122LNm   12 1725    0   0  765 1214 1942 889  680 861

ns
  N01   N02   N03   N04   N05   N06   N07   N08   N09   N10 
 3441 32538  2389   213  3509  4650 23143  5701 13355  2894 

spatial

011LNd

012LNd

020LNm

031LNd

041SOd

042SOd

051LNn

052LNn

061LNd

062LNd

071NPd

072NPd

080LNn

091EYn

092EYn

100LNd

111LNn

112LNn

121LNm

122LNm

nn hm

cd <- colData(sce)[c("lv2")]
df <- data.frame(cd, fq)
mu <- df |>
    group_by(lv2) |>
    summarise(
        across(where(is.numeric),
        \(.) mean(., na.rm=TRUE)))
# scaling
mv <- mu |>
    pivot_longer(-lv2) |>
    group_by(name) |> 
    mutate_at("value", .z)
# wrangling
mx <- pivot_wider(mv)
my <- as.matrix(mx[, -1])
rownames(my) <- mx[[1]]
# plotting
lab <- "z-scaled frequency of X in radial neighborhood of Y"
pal <- colorRamp2(c(-2.5, 0, 2.5), c("navy", "ivory", "maroon"))
set.seed(3); hm <- Heatmap(my, pal, name=lab,
    row_title=NULL, 
    column_title=NULL,
    width=unit(10, "cm"), 
    height=unit(10, "cm"),
    km=10, column_km=10,
    row_dend_width=unit(0.5, "cm"),
    column_dend_height=unit(0.5, "cm"),
    heatmap_legend_param=list(
        direction="horizontal",
        grid_height=unit(2, "mm"),
        legend_width=unit(4, "cm"),
        at=. <- seq(-2.5, 2.5, 0.5),
        title_position="topcenter",
        labels=ifelse(. %% 2 == 0, ., "")))
draw(hm, heatmap_legend_side="top")

fq by sid

ps <- list(
    .plt_fq(sce, "ctx", "sid") + scale_fill_manual(values=.pal_sid),
    .plt_fq(sce, "sid", "ctx") + scale_fill_manual(values=.pal_ctx))
wrap_plots(ps) & coord_equal(20, expand=FALSE)

fq by lv1

ps <- list(
    .plt_fq(sce, "ctx", "lv1") + scale_fill_manual(values=.pal_sub),
    .plt_fq(sce, "lv1", "ctx") + scale_fill_manual(values=.pal_ctx))
wrap_plots(ps) & coord_equal(20, expand=FALSE)

fq by lv2

ps <- list(
    .plt_fq(sce, "ctx", "lv2") + scale_fill_manual(values=.pal_kid),
    .plt_fq(sce, "lv2", "ctx") + scale_fill_manual(values=.pal_ctx))
wrap_plots(ps) & coord_equal(20, expand=FALSE)

cd <- colData(sce)[c("lv2")]
df <- data.frame(fq, cd)
fd <- pivot_longer(df, -names(cd))
mu <- fd |>
    group_by(lv2, name) |>
    mutate_at("name", factor) |>
    summarize_at("value", mean)
mv <- filter(mu, name != "tum")
xo <- mu |>
    filter(name == "tum") |>
    arrange(value) |> 
    pull(lv2)
ggplot(mu, aes(lv2, value, fill=name)) + 
    ggtitle("cells in radial neighborhood of x") +
ggplot(mv, aes(lv2, value, fill=name)) + 
    ggtitle("excluding tumor fraction") +
plot_layout(guides="collect") &
    geom_col(position="fill", key_glyph="point", show.legend=TRUE) &
    scale_fill_manual(NULL, drop=FALSE, values=.pal_kid) &
    scale_x_discrete(limits=xo) &
    .thm_fig_d("minimal", "f") & 
    coord_equal(expand=FALSE) &
    theme(aspect.ratio=2/3,
        axis.title=element_blank(),
        axis.text.y=element_blank(),
        axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))

appendix

save data

saveRDS(fq, "outs/fqs.rds")
saveRDS(ns, "outs/ctx.rds")

session

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Madrid
tzcode source: internal

attached base packages:
[1] stats4    grid      stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] patchwork_1.3.2             ggplot2_4.0.0              
 [3] SingleCellExperiment_1.30.1 SummarizedExperiment_1.38.1
 [5] Biobase_2.68.0              GenomicRanges_1.60.0       
 [7] GenomeInfoDb_1.44.3         IRanges_2.42.0             
 [9] S4Vectors_0.46.0            BiocGenerics_0.54.1        
[11] generics_0.1.4              MatrixGenerics_1.20.0      
[13] matrixStats_1.5.0           ComplexHeatmap_2.24.1      
[15] circlize_0.4.16             cluster_2.1.8.1            
[17] tidyr_1.3.1                 dplyr_1.1.4                
[19] arrow_22.0.0                RANN_2.6.2                 

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1        vipor_0.4.7             farver_2.1.2           
 [4] S7_0.2.0                fastmap_1.2.0           digest_0.6.37          
 [7] lifecycle_1.0.4         Cairo_1.7-0             magrittr_2.0.4         
[10] compiler_4.5.1          rlang_1.1.6             tools_4.5.1            
[13] yaml_2.3.10             knitr_1.50              labeling_0.4.3         
[16] S4Arrays_1.8.1          htmlwidgets_1.6.4       bit_4.6.0              
[19] DelayedArray_0.34.1     mapproj_1.2.12          RColorBrewer_1.1-3     
[22] abind_1.4-8             withr_3.0.2             purrr_1.1.0            
[25] colorspace_2.1-2        scales_1.4.0            iterators_1.0.14       
[28] pals_1.10               dichromat_2.0-0.1       cli_3.6.5              
[31] rmarkdown_2.30          crayon_1.5.3            rstudioapi_0.17.1      
[34] httr_1.4.7              rjson_0.2.23            ggbeeswarm_0.7.2       
[37] maps_3.4.3              assertthat_0.2.1        parallel_4.5.1         
[40] ggrastr_1.0.2           XVector_0.48.0          vctrs_0.6.5            
[43] Matrix_1.7-4            jsonlite_2.0.0          GetoptLong_1.0.5       
[46] bit64_4.6.0-1           clue_0.3-66             beeswarm_0.4.0         
[49] magick_2.9.0            foreach_1.5.2           glue_1.8.0             
[52] codetools_0.2-20        shape_1.4.6.1           gtable_0.3.6           
[55] UCSC.utils_1.4.0        tibble_3.3.0            pillar_1.11.1          
[58] htmltools_0.5.8.1       GenomeInfoDbData_1.2.14 R6_2.6.1               
[61] doParallel_1.0.17       evaluate_1.0.5          lattice_0.22-7         
[64] png_0.1-8               Rcpp_1.1.0              SparseArray_1.8.1      
[67] xfun_0.54               pkgconfig_2.0.3         GlobalOptions_0.1.2